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ABSTRACT 

The  singularly  constrained  generalized  network  problem  represents 
a largo  class  of  capacitated  linear-  programming  (LP)  problems.  This 
class  includes  any  LP  problem  whose  coefficient  matrix,  ignoring  single 
upper  bound  constraints,  contains  m + 1 rows  which  may  be  o-derod  such 
that  each  column  has  at  most  two  non-zero  entries  in  the  first  m rows. 

The  paper  describes  efficient  procedures  for  solving  such  problems  and 
presents  computational  results  which  indicate  that  procedures  are  at  least 
five  times  faster  than  the  state  of  the  art  LP  systems  M PS -360  and  APEX- III. 
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1.0  INTRODUCTION 
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The  singularly  constrained  generalized  network  problem  represents 
a large  class  of  capacitated  linear  programming  (LP)  problems.  This  class 
includes  any  LP  problem  whose  coefficient  matrix,  ignoring  simple  upper 
bound  constraints,  contains  m + 1 rows  which  may  be  ordered  such  that 
each  column  has  at  most  two  non-zero  entries  in  the  first  m rows.  There 
are  numerous  real  world  applications  of  problems  containing  this  structure. 
For  example,  a ma.jor  automobile  manufacturer  has  devised  a singularly 
constrained  generalized  network  model  in  which  the  extra  constraint  is  due 
to  tariff  regulations  that  dollar  flows  across  the  U.S.  -Canadian  border  be 
equal.  Another  example  is  a cash  flow  model  [4]  used  by  the  U.S.  Treasury 
Department  to  analyze  the  effect  of  tax  regulations  on  multinational  firms 
in  which  the  extra  constraint  is  used  to  balance  direct  loans. 

A large  portion  of  the  literature  on  LP  problems  has  been  devoted 
to  special  cases  of  the  singularly  constrained  generalized  network  problem. 
In  particular,  if  the  m + 1st  row  of  the  coefficient  matrix  is  a zero  vector, 
then  depending  on  the  structure  of  the  first  m rows,  the  problem  is  either 
a shortest  path,  assignment,  transportation,  transshipment,  generalized 
transportation,  or  generalized  transshipment  problem.  The  most  efficient 
procedures  for  solving  these  specializations  are  based  on  viewing  the 
problems  in  a graphical  context.  In  particular,  the  simplex  algorithm  has 
been  adapted  [1,  2,  7,  10,  14,  15,  18,  20,  21,  22,  25]  to  solve  problems 
in  which  the  coefficient  matrix  and  the  basis  matrix  are  stored  as  graphs 


using  computer  list  structures. 


The  use  of  such  structures  reduces  both  the  amount  of  work  needed 
to  perform  the  basic  simplex  steps  and  the  amount  of  computer  memory 
required  to  store  essential  data.  The  matrix  operations  of  finding  the 
representation  of  an  entering  vector  and  determining  updated  pseudo  dual 
variable  values  are  performed  by  tracing  paths  within  the  basis  graph. 

Since  the  graphs  contain  only  the  non-zero  entries  in  the  problem  (basis), 
the  list  procedures  [15,  23,  24]  eliminate  checking  or  performing  unnecessary 
arithmetic  operations  on  zero  elements.  In  addition,  by  exploiting  the 
triangular  or  near-triangular  properties  of  such  problems,  these  procedures 
allow  the  basis  matrix  inverse  to  be  stored  implicitly  as  a basis  graph. 

This  graph  is  updated  during  basis  exchange  steps  by  simply  changing  a 
few  pointers  in  the  list  structures.  This  entirely  eliminates  the  arithmetic 
operations  normally  used  to  update  the  basis  inverse  and  also  eliminates 
possible  round-off  error. 

The  advantages  of  such  procedures  are  dramatically  illustrated  by 
the  fact  that  such  simplex  codes  [8]  for  solving  capacitated  transshipment 
problems  which  use  graphical  processing  techniques  [7]  are  at  least  75 
times  faster  than  state  of  the  art  commercial  LP  codes.  In  fact,  a graphical 
processing  code  for  solving  transportation  problems  [9]  has  been  shown  to 
be  150  times  faster  than  the  OPHELIE  LP  code. 

In  addition  to  improving  solution  speed,  the  network  processing 
techniques  have  the  added  advantage  of  requiring  less  computer  memory  to 
solve  a problem.  This  allows  larger  problems  to  be  solved  without  resorting 
to  slower  external  storage  procedures.  If  the  use  of  external  storage  is 
unavoidable,  the  graphical  procedures  are  normally  able  to  keep  the  entire' 
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basis  in  central  memory  and  to  facilitate  efficient  input/output  processing 
of  non-basic  variable  data. 

Motivated  by  these  advantages,  this  paper  develops  efficient  procedures 
for  solving  singularly  constrained  generalized  network  problems  via  the 
simplex  method.  These  procedures  exploit  the  topological  structure  of  the 
first  m rows  of  the  coefficient  matrix.  In  particular,  the  working  basis  is 
partitioned  so  that  the  near  triangularity  property  of  the  basis  can  be  fully 
exploited  from  both  the  algebraic  and  graphical  processing  viewpoints. 

The  importance  of  this  research  is  actually  twofold.  First,  it  has 
resulted  in  a highly  efficient  algorithm  and  computer  code.  Computational 
testing  of  the  proposed  algorithm  indicates  that  it  is  five  times  more  efficient 
than  the  general  purpose  LP  codes,  MPS-360  and  APEX-111.  Second,  it 
focuses  attention  of  the  importance  and  practicality  of  viewing  many  LP 
problems  as  being  network  related,  especially  during  the  formulation  stage's. 


2. 


0 PROBLEM  DEFINITION 

A singularly  constrained  generalized  network  problem  is  defined  as  follows: 
T T 

Minimize  Cjx  + c^y  (1) 

subject  to: 


Gx  + Oy  = b (2) 

fix  + hy  f0  (3) 

0 s x £ u (4) 

0 s y S V (5) 

v/here  each  column  of  G contains  at  most  two  non-zero  entries.  The 


constraints  (2)  will  henceforth  be  referred  to  as  the  underlying  generalized 
network  problem  and  equation  (3)  as  the  extra  constraint. 


1 
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It  is  often  useful  in  the  process  of  formulating  and  solving  uneonsi  raineil 
generalized  network  problems  to  view  them  in  a general  graph  structure  16 j. 

A general  graph  is  an  extension  of  a simple  graph  G(V,  K),  where  V is  a 
finite  non-empty  set  of  vertices  and  E is  a finite  set  of  edges,  each  of  which 
is  identified  with  an  unordered  pair  of  distinct  elements  of  V.  Each  edge  of 
E connects  two  vertices  which  are  then  considered  adjacent.  If  E is  expanded 
fo  contain  edges  that  have  both  endpoints  incident  on  the  same  vertex  (loops) 
or  multiple  edges  connecting  the  same  two  vertices  (parallel  edges)  then 
G(V,E)  is  called  a general  graph. 

The  underlying  generalized  network  of  the  problem  (1-5)  forms  such 
a general  graph.  Each  column  of  G having  two  non-zero  entries  relates  to 
an  edge  connecting  two  vertices.  Each  column  containing  a single  non-zero 
entry  forms  an  edge  with  both  endpoints  incident  on  the  same  vertex.  Such 
edges  will  he  referred  to  as  loops.  Associated  with  each  edge  is  a variable, 
an  upper  bound,  a cost  coefficient,  and  two  non-zero  mulfipliers  (one  in  the 
case  of  loops).  The  multipliers  govern  how  much  of  the  variable  is  to  be 
added  to  or  subtracted  front  the  appropriate  vertex  values  (right-hand  sid. 
values). 

Most  of  the  literature  or,  generalized  network  problems  assumes  ! hat 
each  edge  has  exactly  two  non-zero  multipliers  which  arc  of  opposing  sign. 

Further,  it  is  assumed  without  loss  of  generality  that  one  of  the  multipliers 
has  a value  of  -1.  For  this  special  class  of  generalized  network  problems  a 
direction  may  be  imputed  to  each  edge.  The  graph  is  then  called  a digraph  and 
the  variables  are  customarily  considered  to  be  flows  across  the  directed 
edges  or  arcs. 
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If  artificial  variables  are  appended  to  the  problem,  it  may  be  assumed 
without  loss  of  generality  that  the  underlying  generalized  network  problem 
contains  a coefficient  matrix  which  has  full  row  rank.  If  ihe  coefficient  matrix 
does  not  have  full  row  rank,  it  is  possible  to  further  simplify  the  operations 
presented  subsequently  based  on  the  fact  that  any  generalized  network  without 
full  row  rank  is  equivalent  to  a disjoint  set  of  minimum  cost  flow  networks  fll]. 


3.0  BASIS  STRUCTURE 


Let  P.  = 


P. 

1 

fi 


Let  m be  the  number  of  rows  in  G.  Using  the  standard  bounded 
variable  simplex  algorithm,  a basis  for  the  singularly  constrained  generalized 
network  problem  will  be  a linearly  independent  set  of  m + 1 column  vectors 
taken  from  the  coefficient  matrix  corresponding  to  constraints  (2)  and  (4). 

, i = 1,  2,  . . . , m,  m+1  denote  these  basis  vectors,  where  P. 

represents  the  first  m components  and  fj  the  last  component  of  P..  It  is 
important  to  observe  that  the  rank  of  S = {Pj,  P2,  ■ • • . Pm>  Pm+i  ^ is 
exactly  m.  For  notational  convenience,  assume  the  set  consisting  of  the 
first  m vectors  of  S in  linearly  independent.  Letting  Bn  = [ Pj,  P2,  • • • . P1T1) 
and  F = [ f2,  ....  fml  , any  basis  for  the  problem  (1-5)  may  be  partitioned 

as  follows: 


B 


m + 1 


rm  + l 


(6) 


Since  Bn  is  non-singular  and  is  composed  of  vectors  taken  from  (2),  it 
forms  a working  basis  for  the  underlying  generalized  network  problem.  The 
graphical  structure  of  Bn  has  been  characterized  in  the  literature  [5.  13,21,221. 
It  contains  all  of  the  vertices  but  only  a subset  of  the  edges  of  the  original 
general  graph,  and,  hence,  is  called  a spanning  subgraph.  Due  to  the 
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non-singularity  of  Bn,  the  spanning  subgraph  is  simply  a finite  set  of  disjoint 

quasi-trees.  Each  quasi-tree  is  a simple  tree  to  which  a single  edge  has  been 

I 

added.  The  additional  edge  creates  exactly  one  cycle  (a  circular  path)  in 
the  quasi-tree.  It  should  be  noted  that  a loop,  as  defined  previously,  is  a 
cycle  involving  a single  vertex  and  a single  edge. 

An  efficient  method  for  storing  the  general  graph  B has  been  developed, 
called  the  extended  augmented  predecessor  index  (EAPI)  method  [13|.  The 
EAPI  method,  based  on  the  triple-label  method  of  Johnson  [16,  17],  emph  ys 
computer  list  structures  to  represent  the  vertices  and  edges  of  Bn.  These 
lists  provide  the  means  for  performing  the  two  types  of  quasi-tree  traversal 
needed  to  execute  the  basic  simplex  steps.  Each  quasi-tree  may  be  conceptually 
oriented  with  the  cycle  at  the  top  and  all  attached  subtrees  hanging  down  from 
this  "rooted  cycle".  Representations  of  entering  edges  are  found  by  tracing 
paths  from  given  vertices  up  to  and  around  the  associated  cycles.  Pseudo  dual 
variable  values  are  updated  by  locating  all  vertices  situated  "below"  a given 
vertex.  The  efficiency  of  the  EAPI  method  derives  from  the  facts  that  only 
non-zero  basis  coefficients  are  stored,  that  operations  involving  a basis 
matrix  inverse  may  be  performed  using  the  basis  graph,  and  that  basis 
exchange  steps  involve  only  the  updating  of  list  structure  pointers. 

The  simplex  procedures  in  the  following  sections  are  based  on  the 
basis  partitioning  of  (6).  The  efficiency  of  these  procedures  can  best  be 
realized  if  B is  stored  via  the  EAPI  method.  Therefore,  each  operation  will 
correspondingly  be  described  in  terms  of  the  two  quasi-tree  traversal  methods 
briefly  described  earlier.  However,  a thorough  knowledge  of  the  EAPI  method 
will  not  be  required  to  understand  the'  mathematical  operations  involved  in 
performing  the  basic  simplex  steps. 


I 

! 


W 


-7- 


4.0  OPTIMALITY  TESTING 

Once  a basis  has  been  selected  for  any  linear  programming  problem,  the 
simplex  method  dictates  that  it  be  checked  for  optimality  by  computing  updated 
objective  function  coefficients.  For  network  related  problems,  as  with  general 
linear  programming  problems,  an  efficient  method  for  performing  this  operation 
is  to  use  complementary  slackness  to  calculate  pseudo  dual  variable  values 
associated  with  the  current  primal  basis  and  determine  if  this  dual  solution  is 
feasible.  In  other  words,  once  pseudo  dual  variable  values  have  been  determined, 
an  updated  objective  function  coefficient  is  equivalent  to  the  difference  between 
the  left  and  right  hand  sides  of  the  associated  dual  constraint. 

The  dual  problem  associated  with  the  singularly  constrained  generalized 

network  problem  (1-5)  is: 

T v-T 

Maximize  w b + Sfp  - o u 

subject  to: 

wTG  + Sfj  - 

Sf2  - ♦ T 

w',  S unrestricted 

Jr,  t * o . 

Pseudo  dual  variable  values  w and  S associated  with  the  current  primal 
working  basis  B may  be  found  using  the  theorem  of  complementary  slackness, 
which  yields  the  following  system  of  equations: 


wTB„  + SFn  * cj 

(10) 

WTP  , + Sf  = c 

' m+1  m + 1 m + i 

(11) 

. T 

\Jj  v 


T 

* c 1 


S Cr 


(7) 


(8) 

(0) 
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1 


rP  , 

where  (c^  , cm+^)  is  an  appropriately  ordered  vector  of  the  objective  function 

coefficients  associated  with  the  column  vectors  of  B. 

If  the  value  of  S in  (10)  is  known,  w can  be  computed  using 


wTB 


T 


SF 


(12) 


n n n 

Assuming  Bn  is  stored  via  the  EAPI  method,  (12)  may  be  solved  using  an 

extension  of  the  "pricing-out"  procedure  for  generalized  networks.  This  method 
associates  each  variable  in  w with  a vertex  in  the  general  graph  Bn.  A pseudo 
dual  variable  value  associated  with  a single  cycle  vertex  of  a quasi -tree  is 
determined  using  the  "cycle  factor"  [12,  13].  The  remaining  pseudo  dual  variable 
values  in  the  quasi-tree  are  determined  by  exploiting  the  near  triangularity  of 
the  system  of  equations  (12)  and  using  the  downward  traversal  capabilities 
of  the  EAPI  method.  Note  that  if  a basis  exchange  has  just  occurred  and  F has 
not  changed  from  its  previous  value,  the  maximum  number  of  pseudo  dual 
variable  values  that  have  to  be  updated  will  be  those  associated  with  a single 
quasi-tree.  This  particular  instance,  which  is  standard  in  generalized  networks, 
is  due  to  the  deletion  and  insertion  of  a single  equation  in  the  disjoint,  near 
triangular  system  (12). 

The  solution  of  (12)  is  dependent  upon  knowing  the  current  value  of  S. 

Its  value  may  be  determined  by  combining  (11)  and  (12)  above  and  performing 
the  following  simple  algebraic  manipulation: 

fcnE'n1  - SFnB;‘  >Pm+1  + Sfm  + , = cm+1 


S(FnB _ P 


-1 

n m+1 


f ) 
‘m+r 


"Xlpm+1 


(cTb-'p 


mil  " cm+l^(InBn  pm+l 


m + i 


f , .). 
m+l 


(13) 


■MMh 


% 

-<)- 

The  vector  Rn*Pni  , used  twice  in  (13),  is  simply  the  representation  of 
Pm  + i in  terms  of  the  current  basis  B for  the  underlying  generalized  network 

problem.  In  other  words,  it  is  a vector  whose  components,  ^j,  i = 1,  2 m, 

satisfy: 


U*ipi 


•V2  + 


+ ,4.  p 

'm  m 


Using  this  fact,  equation  (13)  may  be  restated  as 
m m 

S = (iil<J>iCi  " Cm^l)/(,?1  '\fi  " fm  + l 


(14) 


(15) 


Computationally,  S may  be  found  by  tracing  paths  in  the  general  graph 
B . P is  related  to  an  edge,  not  in  B„,  that  is  incident  on  at  most  two 
vertices.  Using  the  EAPT  method,  the  unique  paths  from  these  vorlices  to 
their  respective  cycle(s)  may  be  found.  These  paths  in  conjunction  with  the 
non-repeating  path(s)  around  the  associated  cycle(s)  generate  only  the  non-zero 


values  of  the  rAj,  i = 1,  2,  ....  m (See  [13]  ).  As  these  values  are  determined. 


the  terms 


and 


P 


m 

£ 

i 1 


^ici 


f 

m+i 


( m-i  ] 


(16) 

(17) 


may  be  easily  computed.  Upon  completion,  S is  determined  as  S - p/0. 


Note  that  if  P^.  , is  associated  with  a y,  variable  (i.e.  . P , . is  the  zero 
m + 1 J k m+1 

vector),  each  eA.  = 0,  i = 1,  2,  . ..,  m.  No  graph  traversal  is  needed  in  this 
case  and  S = cm+i/fm+l* 


Once  pseudo  dual  variable  valuts  have  been  determined,  updated 


objective  function  coefficients  may  be  calculated  as  follows: 
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ci 


T SI 


[ w\  f 


P; 


f; 


- Ci 


v i 


(18) 


or 


c.  = wTPi  + Sf:  - c. 

l 1 1 


v 1 


(1‘)) 


Dual,  and  hence  primal,  optimality  is  achieved  if  the  following  conditions 
hold  for  all  non  basic  variables: 
c.  * 0,  xt  = 0 , or  yi  r 0 


c.  = 0,  x-  o basic 

i l J \ 


Cj  5 o,  Xj  = Uj  or  Vj  = Vj, 


Vi  . 


(20) 


5.0  FINDING  THK  t»E PRESENTATION  OF  A VECTOR  TO  ENTER  THK  BASIS 
If  the  current  basis  is  non-optimal  (i.  e.  , at  least  one  of  the  conditions 
(20)  is  violated),  then  a pivot  eligible  vector  is  selected  and  a basis  exchange 
step  is  performed.  To  perform  this  basis  exchange  requires  the  determination 


of  the  representation  Z 


of  the  entering  vectot 


■f'id 


in  terms  of 


the  current  basis  B.  This  representation  may  be  found  by  solving  the  following 
system  of  equations: 

Bnzr  + pm+lzr  ' pr 


(211 


^n^1”  + fm+izr 


f . 
r 


(22) 


Employing  what  is  essentially  the  double-reverse  method  of  Charnes 
and  Cooper  [3],  (21'  may  be  written  as 

(23) 


Z. 


B"  1 P 
n r 


R^P  z 
n m + 1 r 
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Substituting  (23)  into  (22)  and  applying  algebraic  manipulation  yields  the 
following: 


-1  -1 

F (B  P -HP  z ) + f z 
r m n r n m + 1 r'  m+1  r 


^m+1  * FnBnlpm+i>zr  = fr 


FnBm  Pr 


(F  1 P - f )/(F  R^P 
n n r t n n m+1 


f ). 

m+1 


(24) 


-1 


The  vector  Bn  Pp  in  the  numerator  of  (24)  is  the  representation  of 
Pr  in  terms  of  the  basis  B for  the  underlying  generalized  network  problem. 
Letting  6j  , i - 1,  2,  . . . , m denote  the  components  of  this  vector  yields: 


3iPi  + e2p2  + ...  + Pmpm  = pp. 


}2L  2 


(25) 


Noting  that  the  denominator  of  (25)  is  exactly  0 determined  in  (16),  (24) 
may  be  stated  as 


m 


<S  B jf i - fr  )/  e . 


(26) 


Computationally,  the  EA  PI  method  may  be  employed  lo  efficiently 

determine  the  non-zero  0.  , i = 1,2 m of  (25).  This  is  directly 

analogous  to  the  procedure  described  in  section  4 for  finding  the  repre- 
sentation of  Pm  + | in  terms  of  B^.  As  the  Sj  are  determined,  the  numerator 
of  (26)  may  be  successively  summed.  Upon  completion,  the  previously 
obtained  value  of  9 may  be  used  to  determini  zp. 

Using  the  value  of  zp  as  found  in  (26)  and  the  representation  of  P 

and  , , in  terms  of  B„,  (23)  may  be  stated  in  vector  form  as 
m+1  n J 


/P 1 • *lzr 


P2  " *2Zr 


' 7J 


The  efficiency  of  (27)  derives  from  the  fact  that  all  required  information  has 
been  generated  in  previous  operations.  There  are,  however,  certain  cases 
in  which  the  computations  involved  in  determining  Z may  be  further  reduced. 

In  particular,  if  is  associated  with  a y^,  k e K,  then  3^  = 0, 
i = 1,  2,  . . . , m.  Thus  (26)  reduces  to 


zr  = -fr'9 


and  (27)  becomes 


; <ye  i 


\ C*  / 6/ 

' m ' 

On  the  other  hand,  Pm  + ^ may  be  associated  with  a y^.,  k e K,  in  which  case. 


0,  i = 1,  2,  ....  m. 


zr  ' <fr  - Bifi  + l 


I 

I 


-!:•*- 


6.0  MAIMTAI \'IN(.  \ 1 AIITITIUNK b>_  JjA  SIS 

Once  the  re  present  al  ion  of  Ihr  \iTtor  i*nt«*i,i n.ut  1 1 u - basis  lias  been  ionnd, 
a hounded  variable  ratio  teat  may  he  performed  to  determine  the  vector  1«> 
leave  the  basis.  This  involve:  using  the  representation,  the  current  variable 
values,  and  the  upper  hounds.  Computationally  it  is  advantageous  to  couple 
this  step  with  the  operations  performed  in  section  5 since  only  the  non-zero 
entries  of  "Z  need  to  be  checked.  The  basis  exchange  stop  must  not  In- 
performed  indiscriminately,  however,  in  order  to  preserve  the  basis 
partitioning.  The  partitioning  of  the  basis  is  the'  foundation  upon  which  those 
procedures  are  built.  The  destruction  of  t ha t partitioning  would  necessarily 
invalidate  the  given  solution  procedures. 

There  are  three  possible  cases  which  may  occur  depending  upon  tin- 
vector  selected  to  leave'  the  basis. 

Case  1:  The  vector  to  leave  the  basis  is  . In  this  case,  the 

pivot  may  be  performed  directly.  The*  underlying  basis  graph  B will  remain 
unchanged.  Note  that  the  current  8.  , i 1,2,  ... , m values  may  be  used  as 
the  future  i X,  2,  . . . , m values  in  the  computation  of  S and  in  the  next 


pivot  operation. 

Case  2 The  vector  to  leave  the  basis  is  P f P 


mi  1 


and  8 . = 0 . 


If  the  basis  exchange  step  is  performed  immediately,  B would  no  longer 
be  a basis  for  the  underlying  generalized  network  auu  the  basis  partitioning 
would  be  destroyed.  In  other  words,  ' I’j , Pg,  . . . P(  ^ , P j f j , . . . , Pm,  P,.  1 
would  not  be  a linearly  independent  set.  However,  if  must  he  true  that 
{ Plt  P2,  • • • . Pj_j»  Pj  + 1,  • • • . Pm,  Pm  + , 1 is  a linearly  independent  set.  The 


K» 
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partitioning  rmiy,  thus,  be  preserved  by  switching  P.  and  . After 

performing  this  exchange,  the  pivot  may  be  performed  as  in  ease  1. 

Case  3:  The  vector  to  leave  t ho  basis  is  P.  J + j anc^  Pj  ^ 

The  underlying  basis  graph  may  be  modified  by  the  removal  of  P and 
the  insertion  of  Pp.  Note  that  if  ■*  = = 0 in  this  case,  each  , i - 1,  2,  ....  m 
and  S will  remain  unchanged  during  the  basis  exchange.  Since  S does  not 
change,  the  efficient  procedures  presented  in  section  4 may  be  used  to 
update  (as  opposed  to  recalculate)  the  remaining  dual  variable  values. 

In  both  cases  2 and  3 above,  a change  of  pointers  in  the  computer  list 
structures  will  be  necessary.  This  is  fully  described  in  [13]. 


7.0  NUMERICAL  EXAMPLE 

Consider  the  following  singularly  constrained  generalized  network 
problem: 

Minimize  lx I2+6x 13+5x24+Bx25 M*23+4*34 "‘*„3+’*45+6*4G+0y j+My., 
subject  to 

15 

- -10 
0 

= 0 

“"34  ’"43  "43  "40 

•t 

2x, 


lx12f2x13 


lx 


12  iX24 


ix  -4x2f)  + lx23 


lx 


13 


1 

-—x 

4 24 


lx23  * 4x34  + lx43 

- 3x34  + lx43  + lx45  + 3x46 


25 


+ iix 

4 45 


3x46 


lx 


12 


+ lx 


23 


+ 

45 


- !yi  + iv. 


X.  a 0 


y,-  >'v  * 0 


where  y^  is  an  artificial  variable. 
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An  optimal  basis  for  th<  underlying  generalized  network  problem  ,s 
composed  of  Hu-  vectors  1’^.,,  1’.,^,  I * , 1’.^,  and  l’^g.  Adding  to 

these  the  vector  associated  with  y and  expanding  yields  a basis  for  the 
constrained  problem 


2 

0 

0 

0 

0 

0 I 

1 

I 0 \ 

\ 
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1 n 

P 

m i 1 
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0 
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0 

0 

1 
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1 

; 0 

0 

0 

0 

0 

1 

3 
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1 

\ 0 

0 

0 

1 

0 

1 

The  graph  corresponding  to  13  is  shown  in  Figure  1.  The  numbers  at  the 

endpoints  of  each  edge  are  the  corresponding  non-zero  multipliers.  Not 

3 

indicated  on  the  graph  is  y9  — ■ . 

The  first  step  of  the  algorithm  is  to  compute  a value  for  S.  Letting 
P be  the  vector  associated  with  v . (15)  reduces  to  S M. 

Y2  " 2 

The  system  of  equations  associated  with  (12)  is 


2w  j 

lWo  = 

6 

(2b) 

- lw2 

1 

?w4  " 

5 

(30) 

-4w9 

2w- 

a 

8 

(31) 

lw2 

- 1 w = 

4 • 1 M 

(32) 

-4w^ 

- :^w4 

4 

(33) 

3w4 

• 3w„ 

n 

8 

(34) 
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Concentrating  on  the  equations  involved  in  the  cycle  (i.e.  , (30),  (33),  and  (32)) 
the  method  in  f 1 3 ] may  be  used  to  find  w9  (16  - M)/(l-3)  - -9  t ^ M. 

One  pass  through  the  graph  determines  the  other  dual  variable  values  as 
Wj  = 19/2  - 3/4M,  w3  - -13  + 3/2M,  16  - 2M,  wf  -14  + M,  and 

Wfj  = -14  + 2M.  The  dual  constraint  associated  with  the  non-basie  variable 
x,9  is  violated.  Thus  Pj9  is  selected  as  the  entering  vector. 

Using  the  simplifications  discussed  at  the  end  of  section  5 and 
the  procedure  developed  in  [13]  to  find  t lie  representation  of  Pj9  yields 

1 \ 

2 \ 

\ 

\ 

-3  l 
4 1 

i 

0 

zi2  = f ; 

I ! 

16 

0 

V 

3 

' 4 

The  ratio  test  indicates  that  P^9  is  the  vector  to  leave  the  basis.  Thus  a 

case  1 basis  exchange  is  performed.  The  basis  graph  of  Figure  1 will  not 

change  except  for  the  variable  values.  The  new  values  (found  by  a standard 

5 3 

pivot  operation)  are  x^^  = 7,  2,  x9r.  = , x93  ~ 1,  x^  - T , 

5 

x46  " 3 * and  x12  * U 

Since  the  representation  of  Pj9  is  known  in  terms  of  B , the  calculation 

1 4 

of  S is  simplified.  Setting  0 - 1 / < ^-  - 1)  ~ - 3 and  p = 1/2  - 1 = - 1/2 
yields  S = 2/3.  Using  the  same  procedure  as  discussed  earlier,  the  dual 
variable  values  are  found  to  be  Wj  - 9,  w9  - -23/3,  w,^  • -12,  w4  44/3, 

I 

k 


w-  = -40/3,  and 


■38/3.  This  is  a dual  feasible  solution,  and 


therefore,  the  primal  and  dual  solutions  are  optimal. 

8.0  COMPUTATION  At-  RESULTS 
Code  Development 

A computer  code  was  developed  for  solving  singularly  constrained 
generalized  network  problems  in  order  to  confirm  the  relative  efficiency 
of  advanced  specialized  procedures.  The  code  (NETSG)  utilizes  the  extended 
API  method  and  the  procedures  presented  in  the  preceding  sections.  This 
program  is  written  entirely  in  standard  FORTRAN  IV  and  was  initially 
debugged  and  tested  using  the  RIJN  compiler  on  a CDC  6600  computer. 

Under  those  conditions,  it  occupies  a total  of  11N  + 5A  + D100  words  of 
central  memory,  where  N is  the  number  of  nodes  and  A is  the  number  of 
arcs  for  the  specific  problem.  The  code  was  also  subsequently  tested  on 
an  AMDAHL  470  computer  using  the  IBM  FORTRAN  G compiler  and  on  a 
CDC  CYBER  74  computer  using  the  CDC  FTN  compiler. 

Each  of  the  major  functions  of  the  program  were  separated  and 
placed  in  subroutines.  This  made  debugging  easier  and  a’lowed  for  the  use 
of  timing  procedures  for  each  task.  Such  considerations  are  important  if 
alternative  procedures  or  parametric  values  are  to  be  tested.  The  two 
decision  rules  to  be  considered  in  this  paper  are  candidate  list  size  and 
Big-M  values. 

One  of  the  major  considerations  in  any  simplex  based  computer  code 


is  the  choice  of  a vector  to  enter  the  basis  on  each  iteration.  One  of  the 


most  efficient  procedures  for  achieving  this  is  the  use  of  candidate  lists 
(or  multiple  pricing).  Candidate  lists  decrease  the  amount  of  time  required 
to  search  for  the  "best"  pivot  arc  by  providing  a selection  list  that  can  be 


-18- 

aceessed  after  each  iteration.  A more  complete  description  of  an  S - R 
candidate  list  is  as  follows: 

Examine*  sequentially  each  of  the  vertices  for  the  underlying 
generalized  network.  For  each  vertex,  select  the  non  basic  edge 
that  violates  dual  feasibility  the  most  (20).  Add  this  edge  (if  one 
exists)  to  the  list  and  proceed  to  the  next  vertex.  This  is  done  until 
the  list  contains  R entries  or  until  all  of  the  vertices  have  been 
explored.  From  the  list  the  "best  ' edge  is  selected  to  enter  the 
basis.  The  list  is  accessed  after  each  Divot  until  it  is  void  of 
eligible  edges  or  until  it  has  been  accessed  S times.  At  such  time, 
the  list  is  refilled  by  again  examining  each  vertex,  starting  at  i.he 
last  examined  vertex.  When  the  last  (largest  numerical  designation) 
is  encountered,  the  non-graph  dual  constraints  (9)  are  checked  before- 
proceeding  to  the  first  vertex.  If  a complete  pass  through  the  vertices 
is  made  without  filling  the  list,  the  size  of  R is  reset  to  the  number 
found  and  S is  reduced  to  1/2R. 

This  type  of  procedure  has  been  found  t.o  be  exceptionally  efficient  for  pure 
and  generalized  network  codes.  For  this  reason,  the  effect  of  the  initial 
values  of  S and  R for  the  singularly  constrained  code  will  be  shown  subsequently. 

NETSG  currently  implements  an  artificial  start  procedure.  In  other 
words,  the  basis  for  the  underlying  generalized  network  is  composed  entirely 
of  self-loops  with  Rig-M  cost  coefficients.  This  method,  although  probably 
not  the  most  efficient,  provides  a very  quick  initial  basis.  However,  the 
magnitude  of  Big-M  must  be  selected.  This  can  be  a very  critical  determination 


i 


I 


both  in  terms  of  total  solution  time  and  total  number  of  pivots  performed. 
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Choice  of  Decision  Rules 

The  choice  of  the  correct  candidate*  list  size  and  the  magnitude  of 
Big-M  were  both  obtained  by  performing  extensive  computational  tests. 

This  testing  required  the  use  of  a program  to  generate  problem  sets. 
Klingman,  Napier,  and  Stutz  have  created  the  program  NETGEN  for 
generating  pure  network  problems  [19].  With  NETGEN,  the  user  can 
specify  the  desired  structure  of  a problem  and  it  will  be  generated  using 
a random  function.  For  test  pvirposes,  NETGEN  was  modified  to  create 
generalized  arcs  (non-unity  multipliers)  and  coefficients  for  the  extra 
constraint.  Both  of  tv^se  are  created  using  ranges  specified  by  the  user. 

The  input  specifications  for  the  first  three  test  problems  can  be  found  in 
table  1.  All  of  these  problems  have  1000  nodes,  but  thev  differ  dramatically 
in  both  the  structure  of  the  underlying  generalized  network  and  of  the  extra 
constraint. 

The  first  test  was  to  determine  the  proper  initial  size  of  the  candidate 
list.  The  results  of  this  test  can  be  found  in  table  2.  The  sizes  tested  were 
1-1,  5-10,  8-16,  and  15-25.  In  all  three  of  the  problems,  the  8-16  list 
was  clearly  the  best  choice.  On  problem  2,  the  8-16  list  proved  to  be  over 
twice  as  fast  as  the  1-1  list  procedure. 

The  next  decision  rule  to  be  tested  was  the  magnitude  of  Big-M. 

In  the  previous  test,  this  value  had  been  set  to  1000.  Additional  values 
tested  were  500,  250,  150,  and  100.  Each  of  these  were  coupled  with  an 
8-16  candidate  list.  The  results  can  be  seen  in  table  ".  in  problems  1 and 
3,  the  150  value  proved  to  be  the  best  choice.  In  problem  2,  the  250  value 
yielded  the  best  total  time  but  the  150  value  gave  the  lowest  number  of  total 


CANDIDATE  LIST  SIZE 


BIG-M  VALUES 


pivots.  Note  that  the  100  value  (equal  the  largest  cost)  proved  to  yield 
infeasible  solutions  in  all  three  cases. 

Code  Comparisons 

The  next  phase  of  testing  was  to  compare  the  singularly  constrained 
code  NETSG  with  two  widely  used  standard  linear  programming  packages. 

In  this  case,  the  five  problems  indicated  in  table  4 were  used  as  test  data. 
These  problems  range  from  a 50  x 50  constrained  generalized  transportation 
problem  to  a 1000  node  constrained  transshipment  problem. 

The  first  comparison  was  made  against  the  IBM  MPS-360  linear 
programming  package.  The  tests  were  performed  on  an  AMDAHL  470 
computer.  Since  both  IBM-360  and  AMDAHL  470  computers  have  smaller 
word  sizes  than  a CDC  6600,  NETSG  would  not  initially  run  due  to  precision 
errors.  Thus  all  type  real  variables  were  changed  to  double  precision. 

The  effect  of  this  was  two-fold.  First,  it  greatly  increased  the  amount  of 
central  memory  required  to  run  the  program.  Second,  the  use  of  double 
precision  arithmetic  routines  significantly  increased  computation  time. 

The  results  of  testing  against  MPS-360  can  be  found  in  Table  5. 

Even  with  the  disadvantages  listed  above,  NETSG  proved  to  be  between 
5 and  8 times  faster  than  MPS-360.  The  efficiency  of  NETSG  was  much 
more  apparent  as  the  problem  size  increased. 

The  last  comparison  was  performed  using  the  APEX-III  linear 
programming  package.  This  program  is  distributed  by  Control  Data 
Corporation  and  the  tests  were  conducted  on  a CDC  CYBER -74  computer. 

In  this  case,  there  were  no  problems  with  word  size  or  precision.  The 
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basis  of  comparison,  however,  is  not  CPU  seconds,  but  rather  a CDC 
accounting  measure  called  an  SBU.  The  total  SHU  count  for  a job  is 
an  accumulation  of  CPU  seconds  used,  1 10  functions  performed,  and 
central  memory  requirements.  This  measure  seemed  appropriate, 
especially  since  it  allowed  the  comparison  of  actual  dollar  amounts. 

CDC  charges  customers  a minimum  of  $0.  18  per  SBU  used. 

The  results  of  the  comparison  of  NETSG  and  APEX-III  are  contained 
in  Table  6.  In  this  case,  NETSG  was  from  6 to  28  times  more  efficient, 
depending  upon  the  problem  size.  It  should  be  noted  that  in  problem  5 
APEX-III  used  over  20  times  as  many  SBU's  and  was  still  not  at  the 
optimum  after  completing  10,000  iterations. 


i 


NETSG  vs.  M PS -3 60* 


Problem 

NETSG 

MPS-360 

1 

1.  64 

8.  59 

2 

3.  92 

31. 27 

3 

2.  61 

16.48 

4 

13.  06 

110. 68 

5 

35.  62 

1 

DNR 

*times  listed  are  CPI)  seconds 
^problem  data  set  was  too  large 


TABLE  6 


NETSG  vs.  APEX  - III 
NETSG  APEX-III 


Problem 

SBU's1 

T 

Cost 

SBU's 

Cost 

1 

9.  16 

$ 1.65 

62.  70 

$ 11.29 

2 

16.  10 

2.  90 

205. 87 

37.  06 

3 

11.  38 

2.  05 

130. 18 

23.43 

4 

32.  72 

5.  89 

943. 25 

169. 79 

5 

83.  13 

14.  96 

3 

1875. 55 

337.  60 

^System  Billing  Unit 

^Computed  at  $0.  18  per  SBU 
3 

Terminated  after  10,  000  iterations 
Objective  function  value  = 379,065,627 
Optimal  value  = 4,  745,  739 
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